Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(car)
library(tidybayes) # for more tidying outputs
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(rstan) # for interfacing with STAN
library(HDInterval) # for HPD intervals
library(emmeans) # for marginal means etc
library(DHARMa) # for residual diagnostics
library(broom) # for tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
library(patchwork) # for multiple plots
library(ggridges) # for ridge plots
library(bayestestR) # for ROPE
library(see) # for some plots
source("helperFunctions.R")

Scenario

Loyn (1987) modelled the abundance of forest birds with six predictor variables (patch area, distance to nearest patch, distance to nearest larger patch, grazing intensity, altitude and years since the patch had been isolated).

Regent honeyeater

Format of loyn.csv data file

ABUND DIST LDIST AREA GRAZE ALT YR.ISOL
.. .. .. .. .. .. ..
ABUND Abundance of forest birds in patch- response variable
DIST Distance to nearest patch - predictor variable
LDIST Distance to nearest larger patch - predictor variable
AREA Size of the patch - predictor variable
GRAZE Grazing intensity (1 to 5, representing light to heavy) - predictor variable
ALT Altitude - predictor variable
YR.ISOL Number of years since the patch was isolated - predictor variable

The aim of the analysis is to investigate the effects of a range of predictors on the abundance of forest birds.

Read in the data

loyn <- read_csv("../public/data/loyn.csv", trim_ws = TRUE)
glimpse(loyn)
## Rows: 56
## Columns: 7
## $ ABUND   <dbl> 5.3, 2.0, 1.5, 17.1, 13.8, 14.1, 3.8, 2.2, 3.3, 3.0, 27.6, 1.8…
## $ AREA    <dbl> 0.1, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.…
## $ YR.ISOL <dbl> 1968, 1920, 1900, 1966, 1918, 1965, 1955, 1920, 1965, 1900, 19…
## $ DIST    <dbl> 39, 234, 104, 66, 246, 234, 467, 284, 156, 311, 66, 93, 39, 40…
## $ LDIST   <dbl> 39, 234, 311, 66, 246, 285, 467, 1829, 156, 571, 332, 93, 39, …
## $ GRAZE   <dbl> 2, 5, 5, 3, 5, 3, 5, 5, 4, 5, 3, 5, 2, 1, 5, 5, 3, 3, 3, 2, 2,…
## $ ALT     <dbl> 160, 60, 140, 160, 140, 130, 90, 60, 130, 130, 210, 160, 210, …

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ log(\mu_i) = \boldsymbol{\beta} \bf{X_i}\\ \beta_0 \sim{} \mathcal{N}(3,0.5)\\ \beta_{1-9} \sim{} \mathcal{N}(0,2.5)\\ \sigma \sim{} \mathcal{Gamma}(2,1)\\ OR\\ \sigma \sim{} \mathcal{t}(3,0,2.5) \]

where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the additive effects of the scaled versions of distance (ln), distance to the nearest large patch (ln), patch area (ln), grazing intensity, year of isolation and altitude on the abundance of forest birds.

Fit the model

MCMC sampling diagnostics

Model validation

Partial effects plots

Model investigation

ROPE

Region of Practical Equivalence

For standardised parameter (negligible effect)

0.1 * sd(log(loyn$ABUND))
## [1] 0.09024351
loyn.brm3 %>% bayestestR::rope_range()
## [1] -0.01  0.01
loyn.brm3 %>% bayestestR::rope(range = c(-0.09, 0.09))
loyn.brm3 %>%
  bayestestR::rope(range = c(-0.09, 0.09)) %>%
  plot(data = loyn.brm3)

Further analyses

Contrasts

brms

Compare effect of grazing at different patch areas

loyn.list <- with(loyn, list(AREA = c(min(AREA), mean(AREA), max(AREA))))

loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs(reverse = FALSE)
## AREA = 0.1:
##  contrast           ratio lower.HPD upper.HPD
##  fGRAZE1 / fGRAZE2  2.060   0.80527     4.079
##  fGRAZE1 / fGRAZE3  2.254   1.00147     4.150
##  fGRAZE1 / fGRAZE4  9.269   1.14624    48.334
##  fGRAZE1 / fGRAZE5 11.676   0.50891    75.786
##  fGRAZE2 / fGRAZE3  1.087   0.34494     2.247
##  fGRAZE2 / fGRAZE4  4.513   0.40771    25.388
##  fGRAZE2 / fGRAZE5  5.748   0.41880    38.031
##  fGRAZE3 / fGRAZE4  4.066   0.35726    21.807
##  fGRAZE3 / fGRAZE5  5.158   0.25737    35.652
##  fGRAZE4 / fGRAZE5  1.283   0.01567    13.014
## 
## AREA = 69.2696428571429:
##  contrast           ratio lower.HPD upper.HPD
##  fGRAZE1 / fGRAZE2  0.906   0.69666     1.154
##  fGRAZE1 / fGRAZE3  0.849   0.61356     1.089
##  fGRAZE1 / fGRAZE4  0.503   0.22953     0.893
##  fGRAZE1 / fGRAZE5  1.428   0.31054     3.977
##  fGRAZE2 / fGRAZE3  0.941   0.62344     1.261
##  fGRAZE2 / fGRAZE4  0.553   0.22720     1.003
##  fGRAZE2 / fGRAZE5  1.584   0.37070     4.530
##  fGRAZE3 / fGRAZE4  0.591   0.22568     1.056
##  fGRAZE3 / fGRAZE5  1.694   0.37319     4.855
##  fGRAZE4 / fGRAZE5  2.963   0.53841     9.935
## 
## AREA = 1771:
##  contrast           ratio lower.HPD upper.HPD
##  fGRAZE1 / fGRAZE2  0.604   0.29743     1.002
##  fGRAZE1 / fGRAZE3  0.524   0.25328     0.940
##  fGRAZE1 / fGRAZE4  0.117   0.00237     0.487
##  fGRAZE1 / fGRAZE5  0.486   0.00733     3.937
##  fGRAZE2 / fGRAZE3  0.876   0.29757     1.822
##  fGRAZE2 / fGRAZE4  0.196   0.00604     0.846
##  fGRAZE2 / fGRAZE5  0.827   0.01300     6.602
##  fGRAZE3 / fGRAZE4  0.227   0.00640     0.992
##  fGRAZE3 / fGRAZE5  0.961   0.01125     7.687
##  fGRAZE4 / fGRAZE5  4.463   0.02079    62.340
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95
newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  as.data.frame()
head(newdata)
newdata <- loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  pairs() %>%
  gather_emmeans_draws()

newdata %>%
  median_hdci() %>%
  ggplot() +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_pointrange(aes(y = .value, ymin = .lower, ymax = .upper, x = contrast)) +
  facet_wrap(~AREA) +
  coord_flip()

loyn.brm4b %>%
  emmeans(~ fGRAZE | AREA, at = loyn.list, type = "response") %>%
  gather_emmeans_draws()
newdata.p <- newdata %>% summarise(P = sum(.value > 1) / n())
g <- newdata %>%
  ggplot() +
  geom_vline(xintercept = 1, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = contrast,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
  facet_grid(~ round(AREA, 1))
g + geom_text(data = newdata.p, aes(y = contrast, x = 1, label = paste("P = ", round(P, 3))), hjust = -0.2, position = position_nudge(y = 0.5))

Summary figure

References

Loyn, R. H. 1987. “Nature Conservation: The Role of Remnants of Native Vegetation.” In, edited by D. A. Saunders, G. W. Arnold, A. A. Burbridge, and A. J. M. Hopkins. Chipping Norton, NSW: Surrey Beatty & Sons.